sound_field_analysis example 1 - Ideal plane wave simulation

In the following notebook, the sound field of an ideal plane wave impinging on a spherical microphone array is generated and visualized. Then, the impulse response of two different looking direction is extraced using plane wave decomposition.

Setup

In [1]:
import sys
sys.path.insert(0, '../')
from sound_field_analysis import gen, plot, utils, process

from plotly.offline import init_notebook_mode
init_notebook_mode()

Configuration

First the limit for the spherical domain is set. Lower order correspond to faster calculation time at a loss of fidelity in the radial spatial domain.

Then, the configuration of the microphone array that is to be simulated. It simply holds three elements:

  • radius of the array
  • type of the sphere (open, rigid or dual)
  • type of transducers (omni or cardioid)

Furthermore, the wavetype is set to 'plane' and the direction of the wave is specified to come from "straight ahead". NFFT denotes the number of the FFT bins and fs the sampling frequency used for the simulation.

In [2]:
order = 5

array_configuration = [0.5, 'rigid', 'cardioid']

wavetype = 'plane'
azimuth = utils.deg2rad(0)
colatitude = utils.deg2rad(90)

NFFT = 128
fs = 48000

Generation and radial filters

gen.ideal_wave() returns the spherical fourier coefficiens of a ideal plane wave as captured by an ideal array of the provided configuration. To further work with these, it is often helpful to also generate the corresponding radial filters along side.

In [3]:
field_coefficients = gen.ideal_wave(order, fs, azimuth, colatitude, array_configuration, NFFT=NFFT)
radial_filter = gen.radial_filter_fullspec(order, NFFT, fs, array_configuration)

Visualization

To visualize the field, the spherical fourier coefficients and the generated radial filters are passed to plot.makeMTX() along with a kr bin. The function returns the sound pressure at a resolution of 1 degree, which then can be visualized using plot.plot3D(). Several styles ('shape', 'sphere', 'flat') exist.

In [4]:
kr_IDX = 64
vizMTX = plot.makeMTX(field_coefficients, radial_filter, kr_IDX)
plot.plot3D(vizMTX, style='shape')

Extracting impulse responses

In order to extract the impulse response, a plane wave decomposition (PWD) is performed using he process.plane_wave_decomp() function, whos output in the frequency domain can then simply be transformed to the time domain using process.iFFT(). We can calculate several looking directions at once by supplying a vector of azimuth / colatitude pairs.

Since this is an approximated ideal plane wave, the impulse response will be almost dirac impulse-like in the direction of the plane wave (azimuth = 0°, colatitude = 90°) and almost zero at a 90° angle (like azimuth = 90°, colatitude = 90°).

In [5]:
looking_direction = [[0, utils.deg2rad(90)],
                     [utils.deg2rad(90), utils.deg2rad(90)]]
Y = process.plane_wave_decomp(order, looking_direction, field_coefficients, radial_filter)

impulseResponses = process.iFFT(Y)

Plot time signal and frequency response

In [6]:
spectrum, f = process.FFT([impulseResponses, fs])
plot.plot2D(impulseResponses, type='time', fs=fs)
plot.plot2D(utils.db(spectrum), type='logFFT', fs=fs)